Immune population analysis

General Setup

Setup chunk

Setup reticulate

knitr::opts_chunk$set(fig.width = 8)
knitr::opts_knit$set(root.dir = normalizePath(".."))
knitr::opts_knit$get("root.dir")

[1] "/nas/groups/treutlein/USERS/tomasgomes/projects/liver_regen"

Load libraries

library(reticulate)
knitr::knit_engines$set(python = reticulate::eng_python)
py_available(initialize = FALSE)

[1] FALSE

use_python(Sys.which("python"))
py_config()

python:         /home/tpires/bin/miniconda3/bin/python
libpython:      /home/tpires/bin/miniconda3/lib/libpython3.8.so
pythonhome:     /home/tpires/bin/miniconda3:/home/tpires/bin/miniconda3
version:        3.8.3 (default, May 19 2020, 18:47:26)  [GCC 7.3.0]
numpy:          /home/tpires/bin/miniconda3/lib/python3.8/site-packages/numpy
numpy_version:  1.22.1

NOTE: Python version was forced by RETICULATE_PYTHON

Load data (from all cells)

library(Seurat)

Attaching SeuratObject

library(ggplot2)
library(ggrepel)
library(destiny)
library(plyr)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:plyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(RColorBrewer)

Other functions

# Plot correlations
plotCorr = function(cort, sp1, sp2){
  # cluster and order labels
  hcr = hclust(dist(cort$r), method = "ward.D2")
  hcc = hclust(dist(t(cort$r)), method = "ward.D2")
  hcr = hcr$labels[hcr$order]
  hcc = hcc$labels[hcc$order]

  # reshaping the correlations
  plot_df = reshape2::melt(cort$r)
  plot_df$Var1 = factor(plot_df$Var1, levels = rev(hcr))
  plot_df$Var2 = factor(plot_df$Var2, levels = hcc)

  # add pvalue and max cor infor
  plot_df$padj = -log10(reshape2::melt(cort$p.adj+min(cort$p.adj[cort$p.adj>0])/10)$value)
  plot_df$rowmax = apply(Reduce(cbind, lapply(names(cort$maxrow),
                                              function(n) plot_df$Var1==n &
                                                plot_df$Var2==colnames(cort$r)[cort$maxrow[n]])),
                         1, any)
  plot_df$colmax = apply(Reduce(cbind, lapply(names(cort$maxcol),
                                              function(n) plot_df$Var2==n &
                                                plot_df$Var1==rownames(cort$r)[cort$maxcol[n]])),
                         1, any)
  plot_df$markcol = plot_df$value>quantile(plot_df$value, 0.98)

  # getting a colourscale where 0 is white in the middle, and intensity leveled by max(abs(value))
  cols = colorRampPalette(c(rev(RColorBrewer::brewer.pal(9, "Blues")),
                            RColorBrewer::brewer.pal(9, "Reds")))(101)
  br = seq(-max(abs(cort$r)), max(abs(cort$r)), length.out = 101)
  cols = cols[!(br>max(cort$r) | br<min(cort$r))]

  corplot = ggplot()+
    geom_point(data = plot_df, mapping = aes(x = Var2, y = Var1, fill = value, size = padj),
               shape = 21)+
    geom_point(data = plot_df[plot_df$rowmax,], mapping = aes(x = Var2, y = Var1, size = padj),
               shape = "—", show.legend = F, colour = "grey10")+
    geom_point(data = plot_df[plot_df$colmax,], mapping = aes(x = Var2, y = Var1, size = padj),
               shape = "|", show.legend = F, colour = "grey10")+
    scale_x_discrete(expand = c(0,0.7))+
    scale_y_discrete(expand = c(0,0.7))+
    scale_fill_gradientn(breaks = signif(c(min(cort$r)+0.005, 0, max(cort$r)-0.005),2),
                         values = scales::rescale(c(min(br), 0, max(br))),
                         colours = cols)+
    labs(x = sp2, y = sp1, fill = "Spearman's\nrho", size = "-log10\nadj. p-value")+
    theme_classic()+
    theme(axis.title = element_text(colour = "black", face = "bold"),
          axis.text = element_text(colour = "black"),
          axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          legend.title = element_text(size = 9),
          legend.text = element_text(size = 8))

  return(corplot)
}

Immune cell subsetting

Subset immune cells

# Plot correlations
plotCorr = function(cort, sp1, sp2){
  # cluster and order labels
  hcr = hclust(dist(cort$r), method = "ward.D2")
  hcc = hclust(dist(t(cort$r)), method = "ward.D2")
  hcr = hcr$labels[hcr$order]
  hcc = hcc$labels[hcc$order]

  # reshaping the correlations
  plot_df = reshape2::melt(cort$r)
  plot_df$Var1 = factor(plot_df$Var1, levels = rev(hcr))
  plot_df$Var2 = factor(plot_df$Var2, levels = hcc)

  # add pvalue and max cor infor
  plot_df$padj = -log10(reshape2::melt(cort$p.adj+min(cort$p.adj[cort$p.adj>0])/10)$value)
  plot_df$rowmax = apply(Reduce(cbind, lapply(names(cort$maxrow),
                                              function(n) plot_df$Var1==n &
                                                plot_df$Var2==colnames(cort$r)[cort$maxrow[n]])),
                         1, any)
  plot_df$colmax = apply(Reduce(cbind, lapply(names(cort$maxcol),
                                              function(n) plot_df$Var2==n &
                                                plot_df$Var1==rownames(cort$r)[cort$maxcol[n]])),
                         1, any)
  plot_df$markcol = plot_df$value>quantile(plot_df$value, 0.98)

  # getting a colourscale where 0 is white in the middle, and intensity leveled by max(abs(value))
  cols = colorRampPalette(c(rev(RColorBrewer::brewer.pal(9, "Blues")),
                            RColorBrewer::brewer.pal(9, "Reds")))(101)
  br = seq(-max(abs(cort$r)), max(abs(cort$r)), length.out = 101)
  cols = cols[!(br>max(cort$r) | br<min(cort$r))]

  corplot = ggplot()+
    geom_point(data = plot_df, mapping = aes(x = Var2, y = Var1, fill = value, size = padj),
               shape = 21)+
    geom_point(data = plot_df[plot_df$rowmax,], mapping = aes(x = Var2, y = Var1, size = padj),
               shape = "—", show.legend = F, colour = "grey10")+
    geom_point(data = plot_df[plot_df$colmax,], mapping = aes(x = Var2, y = Var1, size = padj),
               shape = "|", show.legend = F, colour = "grey10")+
    scale_x_discrete(expand = c(0,0.7))+
    scale_y_discrete(expand = c(0,0.7))+
    scale_fill_gradientn(breaks = signif(c(min(cort$r)+0.005, 0, max(cort$r)-0.005),2),
                         values = scales::rescale(c(min(br), 0, max(br))),
                         colours = cols)+
    labs(x = sp2, y = sp1, fill = "Spearman's\nrho", size = "-log10\nadj. p-value")+
    theme_classic()+
    theme(axis.title = element_text(colour = "black", face = "bold"),
          axis.text = element_text(colour = "black"),
          axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
          legend.title = element_text(size = 9),
          legend.text = element_text(size = 8))

  return(corplot)
}

Lymphoid cell analysis

Subset Lymphoid cells

immune_pops = c("ab-T cells", "gd-T cells", "Plasmablasts", "B cells")
all_l_cells = allcells_css[,allcells_css@meta.data$allcells_clusters %in% immune_pops]
all_l_cells = suppressWarnings(SCTransform(all_l_cells, do.correct.umi = T, verbose = F, 
                                             vars.to.regress=c("unique_name","nCount_RNA"),
                                             variable.features.rv.th = 1, seed.use = 1,
                                             return.only.var.genes = F, 
                                             variable.features.n = NULL))
all_l_cells = RunPCA(all_l_cells, verbose = F)
all_l_cells = RunUMAP(all_l_cells, dims = 1:25, verbose = F)
DimPlot(all_l_cells, reduction = "umap", group.by = "Condition")
DimPlot(all_l_cells, reduction = "umap", group.by = "allcells_clusters")
DimPlot(all_l_cells, reduction = "umap", group.by = "Donor")
DimPlot(all_l_cells, reduction = "umap", group.by = "Phase")
FeaturePlot(all_l_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "NKG7", "TRGC1", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

Cluster and get markers

all_l_cells = FindNeighbors(all_l_cells, reduction = "pca", dims = 1:25,
                            prune.SNN = 1/5, force.recalc = T, graph.name = "pca25")
all_l_cells = FindClusters(all_l_cells, algorithm = 2, verbose = F, graph.name = "pca25",
                           resolution = seq(0.1, 2, 0.1))
DimPlot(all_l_cells, reduction = "umap", group.by = "pca25_res.0.4", label = T)
DimPlot(all_l_cells, reduction = "umap", group.by = "Donor", label = F)
DimPlot(all_l_cells, reduction = "umap", group.by = "Condition", label = F)

all_l_cells = SetIdent(all_l_cells, value = "pca25_res.0.4")
mk_lcells = FindAllMarkers(all_l_cells, logfc.threshold = 0.2, pseudocount.use = 0.1)
write.csv(mk_lcells[mk_lcells$p_val_adj<=0.05,], 
          file = "results/immune/markers_lymphoid_subpop_all.csv", row.names = T, quote = F)


mk02 = FindMarkers(all_l_cells, ident.1 = "0", ident.2 = "2", 
                   logfc.threshold = 0.2, pseudocount.use = 0.1)

Add annotations

new_l_labs = c("0" = "ab-T cells 1",
               "1" = "NK/gd-T cells",
               "2" = "ab-T cells 2",
               "3" = "Infiltrating NK cells", # PTGDS,CX3CR1 (infilt)
               "4" = "IgA+ Plasma cells",
               "5" = "B cells",
               "6" = "IgG+ Plasma cells",
               "7" = "Dividing NK cells",
               "8" = "ab-T cells (stress)")

all_l_cells$lymphoid_annot = new_l_labs[as.character(all_l_cells$pca25_res.0.4)]

Subset T/NK cells

immune_pops = c("ab-T cells 1", "ab-T cells 2", "NK/gd-T cells", "Infiltrating NK cells")
all_t_cells = all_l_cells[,all_l_cells@meta.data$lymphoid_annot %in% immune_pops]
all_t_cells = suppressWarnings(SCTransform(all_t_cells, do.correct.umi = T, verbose = F, 
                                             vars.to.regress=c("unique_name","nCount_RNA"),
                                             variable.features.rv.th = 1, seed.use = 1,
                                             return.only.var.genes = F, 
                                             variable.features.n = NULL))
all_t_cells = RunPCA(all_t_cells, verbose = F)
all_t_cells = RunUMAP(all_t_cells, dims = 1:25, verbose = F)
DimPlot(all_t_cells, reduction = "umap", group.by = "Condition")
DimPlot(all_t_cells, reduction = "umap", group.by = "allcells_clusters")
DimPlot(all_t_cells, reduction = "umap", group.by = "lymphoid_annot")
DimPlot(all_t_cells, reduction = "umap", group.by = "Donor")
DimPlot(all_t_cells, reduction = "umap", group.by = "Phase")
FeaturePlot(all_t_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "NKG7", "TRGC1", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

Cluster and get markers

all_t_cells = FindNeighbors(all_t_cells, reduction = "pca", dims = 1:25,
                            prune.SNN = 1/5, force.recalc = T, graph.name = "pca25")
all_t_cells = FindClusters(all_t_cells, algorithm = 2, verbose = F, graph.name = "pca25",
                           resolution = seq(0.1, 2, 0.1))
DimPlot(all_t_cells, reduction = "umap", group.by = "pca25_res.0.9", label = T)
DimPlot(all_t_cells, reduction = "umap", group.by = "Donor", label = F)
DimPlot(all_t_cells, reduction = "umap", group.by = "Condition", label = F)
DimPlot(all_t_cells, reduction = "umap", group.by = "t_annot", label = T)

all_t_cells = SetIdent(all_t_cells, value = "pca25_res.0.9")
mk_tcells = FindAllMarkers(all_t_cells, logfc.threshold = 0.2, pseudocount.use = 0.1)
write.csv(mk_tcells[mk_tcells$p_val_adj<=0.05,], 
          file = "results/immune/markers_t_subpop_all.csv", row.names = T, quote = F)

saveRDS(mk_tcells, file = "./results/immune/clust_markers_t.RDS")

mk97 = FindMarkers(all_t_cells, ident.1 = "9", ident.2 = "7", 
                   logfc.threshold = 0.2, pseudocount.use = 0.1)
mk26 = FindMarkers(all_t_cells, ident.1 = "2", ident.2 = "6",
                   logfc.threshold = 0.2, pseudocount.use = 0.1)
mk813 = FindMarkers(all_t_cells, ident.1 = "8", ident.2 = "13",
                   logfc.threshold = 0.2, pseudocount.use = 0.1)

Annotate

# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5007630/
new_t_labs = c("0" = "NK cells 1",
               "1" = "TRM cells", # may have CD4 and CD8; ITGA1
               "2" = "MAIT cells 2", # CXCR6, CCR6, CCR5, some RORC (not DE with the ILC3), some ZBTB16
               "3" = "NK cells 2",
               "4" = "NK cells 3",
               "5" = "CD8 ab-T cells 1",
               "6" = "MAIT cells 1", # hard to be sure, but has many hallmarks, even higher CD8A
               "7" = "Infiltrating NK cells",
               "8" = "Naive CD4+ T cells", 
               "9" = "gd-T cells", # CD3+ vs cl7, TRDC/TRG
               "10" = "CD8 ab-T cells 2",
               "11" = "CD8 ab-T cells (stress)",
               "12" = "Treg",
               "13" = "ILC3") # KIT, AHR, RORC, no CD3/CD8/CD4

all_t_cells$t_annot = new_t_labs[as.character(all_t_cells$pca25_res.0.9)]

Myeloid cell analysis

Subset Myeloid cells

immune_pops = c("cDCs 2", "Macrophages", "Kupffer cells", "cDCs 1", "pDCs")
all_m_cells = allcells_css[,allcells_css@meta.data$allcells_clusters %in% immune_pops]
all_m_cells = suppressWarnings(SCTransform(all_m_cells, do.correct.umi = T, verbose = F, 
                                             vars.to.regress=c("unique_name","nCount_RNA"),
                                             variable.features.rv.th = 1, seed.use = 1,
                                             return.only.var.genes = F, 
                                             variable.features.n = NULL))
all_m_cells = RunPCA(all_m_cells, verbose = F)
all_m_cells = RunUMAP(all_m_cells, dims = 1:25, verbose = F)
DimPlot(all_m_cells, reduction = "umap", group.by = "Condition")
DimPlot(all_m_cells, reduction = "umap", group.by = "allcells_clusters")
DimPlot(all_m_cells, reduction = "umap", group.by = "Donor")
DimPlot(all_m_cells, reduction = "umap", group.by = "Phase")
FeaturePlot(all_m_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "NKG7", "TRGC1", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

Clustering and markers

all_m_cells = FindNeighbors(all_m_cells, reduction = "pca", dims = 1:25,
                            prune.SNN = 1/5, force.recalc = T, graph.name = "pca25")
all_m_cells = FindClusters(all_m_cells, algorithm = 2, verbose = F, graph.name = "pca25",
                           resolution = seq(0.1, 2, 0.1))
DimPlot(all_m_cells, reduction = "umap", group.by = "pca25_res.0.6", label = T)
DimPlot(all_m_cells, reduction = "umap", group.by = "Donor", label = F)
DimPlot(all_m_cells, reduction = "umap", group.by = "Condition", label = F)
DimPlot(all_m_cells, reduction = "umap", group.by = "allcells_clusters", label = F)

all_m_cells = SetIdent(all_m_cells, value = "pca25_res.0.6")
mk_mcells = FindAllMarkers(all_m_cells, logfc.threshold = 0.2, pseudocount.use = 0.1)
write.csv(mk_mcells[mk_mcells$p_val_adj<=0.05,], 
          file = "results/immune/markers_m_subpop_all.csv", row.names = T, quote = F)

saveRDS(mk_mcells, file = "./results/immune/clust_markers_myeloid.RDS")

mk1011 = FindMarkers(all_m_cells, ident.1 = "10", ident.2 = "11",
                     logfc.threshold = 0.2, pseudocount.use = 0.1)

Add annotations

mrks_q = SoupX::quickMarkers(all_m_cells@assays$SCT@counts,
                             all_m_cells@active.ident, N = 10)
View(mrks_q[mrks_q$qval<=0.05,])

new_m_labs = c("0" = "Kupffer cells",
               "1" = "Monocytes/cDCs",
               "2" = "Macrophages",
               "3" = "Monocytes/cDCs",
               "4" = "Kupffer cells",
               "5" = "Monocytes/cDCs",
               "6" = "Monocytes/cDCs",
               "7" = "Macrophages",
               "8" = "Monocytes/cDCs",
               "9" = "cDC1",
               "10" = "pDCs",
               "11" = "pDCs",
               "12" = "Dividing cDCs",
               "13" = "Kupffer cells",
               "14" = "Hepatocytes")

all_m_cells$mye_annot = new_m_labs[as.character(all_m_cells$pca25_res.0.6)]

Subset Monocytes

immune_pops = c("2", "7", "1", "6", "8", "3", "5", "4", "13", "0", "9")
all_mon_cells = all_m_cells[,all_m_cells@meta.data$pca25_res.0.6 %in% immune_pops &
                            all_m_cells@meta.data$allcells_clusters %in% c("Macrophages", "cDCs 1",
                                                                           "cDCs 2", "Kupffer cells")]
all_mon_cells = suppressWarnings(SCTransform(all_mon_cells, do.correct.umi = T, verbose = F, 
                                             vars.to.regress=c("unique_name","nCount_RNA"),
                                             variable.features.rv.th = 1, seed.use = 1,
                                             return.only.var.genes = F, 
                                             variable.features.n = NULL))
all_mon_cells = RunPCA(all_mon_cells, verbose = F)
all_mon_cells = RunUMAP(all_mon_cells, dims = 1:25, verbose = F)
DimPlot(all_mon_cells, reduction = "umap", group.by = "Condition")
DimPlot(all_mon_cells, reduction = "umap", group.by = "allcells_clusters")
DimPlot(all_mon_cells, reduction = "umap", group.by = "pca25_res.0.6")
DimPlot(all_mon_cells, reduction = "umap", group.by = "Donor")
DimPlot(all_mon_cells, reduction = "umap", group.by = "Phase")
FeaturePlot(all_mon_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "NKG7", "TRGC1", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

Clustering and markers

all_mon_cells = FindNeighbors(all_mon_cells, reduction = "pca", dims = 1:25,
                            prune.SNN = 1/5, force.recalc = T, graph.name = "pca25")
all_mon_cells = FindClusters(all_mon_cells, algorithm = 2, verbose = F, graph.name = "pca25",
                           resolution = seq(0.1, 2, 0.1))
DimPlot(all_mon_cells, reduction = "umap", group.by = "pca25_res.0.5", label = T)
DimPlot(all_mon_cells, reduction = "umap", group.by = "allcells_clusters", label = T)
DimPlot(all_mon_cells, reduction = "umap", group.by = "Donor", label = F)
DimPlot(all_mon_cells, reduction = "umap", group.by = "Condition", label = F)

all_mon_cells = SetIdent(all_mon_cells, value = "pca25_res.0.5")
mk_mon = FindAllMarkers(all_mon_cells, logfc.threshold = 0.2, pseudocount.use = 0.1)
write.csv(mk_mon[mk_mon$p_val_adj<=0.05,], 
          file = "results/immune/markers_mon_subpop_all.csv", row.names = T, quote = F)

saveRDS(mk_mon, file = "./results/immune/clust_markers_mon.RDS")

mk61 = FindMarkers(all_mon_cells, ident.1 = "6", ident.2 = "1", 
                   logfc.threshold = 0.2, pseudocount.use = 0.1)
mk21 = FindMarkers(all_mon_cells, ident.1 = "2", ident.2 = "1", 
                   logfc.threshold = 0.2, pseudocount.use = 0.1)

Add annotations

new_m_labs = c("0" = "Kupffer cells",
               "1" = "cDC2",
               "2" = "Monocytes (IGSF21+ GPR34+)", # similar to those identified here https://www.nature.com/articles/s41586-020-2922-4 
               "3" = "Macrophages (HES4+)", # very similar to 5, HES4 is one of its most unique markers
               "4" = "Kupffer cells (SUCNR1+)", # disproves https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1986575/; these KC are inflammatory/antiviral
               "5" = "Macrophages",
               "6" = "Monocytes (TREM2+ CD9+)", # in ramachandran et al, assoc with fibrotic scars
               ## in their projection, also close to KC
               "7" = "cDC2",
               "8" = "cDC1",
               "9" = "Monocytes (secretory)", # related to 2 and 6
               "10" = "activated DCs") # CD80, CD86, CCR7

all_mon_cells$mono_annot = new_m_labs[as.character(all_mon_cells$pca25_res.0.5)]

Put annotations on single immune cell object

Make dataframe with new annotations

newannot_l = list(ldf = all_l_cells@meta.data[,c("allcells_clusters", "lymphoid_annot")],
                  tdf = all_t_cells@meta.data[,c("allcells_clusters", "t_annot")],
                  mdf = all_m_cells@meta.data[,c("allcells_clusters", "mye_annot")],
                  mondf = all_mon_cells@meta.data[,c("allcells_clusters", "mono_annot")])
for(n in names(newannot_l)){
  newannot_l[[n]]$cells = rownames(newannot_l[[n]])
}

newannot_df = Reduce(function(x,y){merge(x,y, by = "cells", all = T)}, newannot_l)[,c(1,3,5,7,9)]

newannot_df$t_annot[is.na(newannot_df$t_annot)] = newannot_df$lymphoid_annot[is.na(newannot_df$t_annot)]
newannot_df$mono_annot[is.na(newannot_df$mono_annot)] = newannot_df$mye_annot[is.na(newannot_df$mono_annot)]
newannot_df$immune_annot = newannot_df$t_annot
newannot_df$immune_annot[is.na(newannot_df$immune_annot)] = newannot_df$mono_annot[is.na(newannot_df$immune_annot)]

newannot_df = data.frame(row.names = newannot_df$cells, 
                         immune_annot = newannot_df$immune_annot)

Add to the immune Seurat object

all_imm_cells = AddMetaData(all_imm_cells, metadata = newannot_df)
# the original "Dividing cells" will be relabeled "Dividing T/NK cells",
## and the newly annotated "Dividing NK cells" will be renamed to match this
all_imm_cells$immune_annot[is.na(all_imm_cells$immune_annot)] = "Dividing T/NK cells"
all_imm_cells$immune_annot[all_imm_cells$immune_annot=="Dividing NK cells"] = "Dividing T/NK cells"
all_imm_cells$immune_annot[all_imm_cells$immune_annot=="Hepatocytes"] = "Hepatocyte-Monocyte interaction"

DimPlot(all_l_cells, group.by = "lymphoid_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Lymphoid cells")
DimPlot(all_t_cells, group.by = "t_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("T and NK cells (subset of Lymphoid)")
DimPlot(all_m_cells, group.by = "mye_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Myeloid cells")
DimPlot(all_mon_cells, group.by = "mono_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Monocytes cells (subset of Myeloid)")
DimPlot(all_imm_cells, group.by = "allcells_clusters", reduction = "umap", label = T)+NoLegend()+ggtitle("Immune cells (original annotation)")
DimPlot(all_imm_cells, group.by = "immune_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Immune cells (new detailed annotation)")

Save Seurat objects

saveRDS(all_l_cells, file = "results/immune/all_l_cells.RDS")
saveRDS(all_t_cells, file = "results/immune/all_t_cells.RDS")
saveRDS(all_m_cells, file = "results/immune/all_m_cells.RDS")
saveRDS(all_mon_cells, file = "results/immune/all_mon_cells.RDS")
saveRDS(all_imm_cells, file = "results/immune/all_imm_cells.RDS")

Analysis of immune cell types

Load data

all_l_cells = readRDS(file = "results/immune/all_l_cells.RDS")
all_t_cells = readRDS(file = "results/immune/all_t_cells.RDS")
all_m_cells = readRDS(file = "results/immune/all_m_cells.RDS")
all_mon_cells = readRDS(file = "results/immune/all_mon_cells.RDS")
all_imm_cells = readRDS(file = "results/immune/all_imm_cells.RDS")

Markers

all_l_cells = readRDS(file = "results/immune/all_l_cells.RDS")
all_t_cells = readRDS(file = "results/immune/all_t_cells.RDS")
all_m_cells = readRDS(file = "results/immune/all_m_cells.RDS")
all_mon_cells = readRDS(file = "results/immune/all_mon_cells.RDS")
all_imm_cells = readRDS(file = "results/immune/all_imm_cells.RDS")

Plot UMAPs

mk_lcells = read.csv("results/immune/markers_lymphoid_subpop_all.csv", 
                     row.names = 1, header = T)
new_l_labs = c("0" = "ab-T cells 1",
               "1" = "NK/gd-T cells",
               "2" = "ab-T cells 2",
               "3" = "Infiltrating NK cells", # PTGDS,CX3CR1 (infilt)
               "4" = "IgA+ Plasma cells",
               "5" = "B cells",
               "6" = "IgG+ Plasma cells",
               "7" = "Dividing NK cells",
               "8" = "ab-T cells (stress)")
mk_lcells$annot = new_l_labs[as.character(mk_lcells$cluster)]

mk_tcells = read.csv("results/immune/markers_t_subpop_all.csv", 
                     row.names = 1, header = T)
new_t_labs = c("0" = "NK cells 1",
               "1" = "TRM cells", # may have CD4 and CD8; ITGA1
               "2" = "MAIT cells 2", # CXCR6, CCR6, CCR5, some RORC (not DE with the ILC3), some ZBTB16
               "3" = "NK cells 2",
               "4" = "NK cells 3",
               "5" = "CD8 ab-T cells 1",
               "6" = "MAIT cells 1", # hard to be sure, but has many hallmarks, even higher CD8A
               "7" = "Infiltrating NK cells",
               "8" = "Naive CD4+ T cells", 
               "9" = "gd-T cells", # CD3+ vs cl7, TRDC/TRG
               "10" = "CD8 ab-T cells 2",
               "11" = "CD8 ab-T cells (stress)",
               "12" = "Treg",
               "13" = "ILC3") # KIT, AHR, RORC, no CD3/CD8/CD4
mk_tcells$annot = new_t_labs[as.character(mk_tcells$cluster)]

mk_mcells = read.csv("results/immune/markers_m_subpop_all.csv", 
                     row.names = 1, header = T)
new_m_labs = c("0" = "Kupffer cells",
               "1" = "Monocytes/cDCs",
               "2" = "Macrophages",
               "3" = "Monocytes/cDCs",
               "4" = "Kupffer cells",
               "5" = "Monocytes/cDCs",
               "6" = "Monocytes/cDCs",
               "7" = "Macrophages",
               "8" = "Monocytes/cDCs",
               "9" = "cDC1",
               "10" = "pDCs",
               "11" = "pDCs",
               "12" = "Dividing cDCs",
               "13" = "Kupffer cells",
               "14" = "Hepatocytes")
mk_mcells$annot = new_m_labs[as.character(mk_mcells$cluster)]

mk_monocells = read.csv("results/immune/markers_mon_subpop_all.csv", 
                        row.names = 1, header = T)
new_m_labs = c("0" = "Kupffer cells",
               "1" = "cDC2",
               "2" = "Monocytes (IGSF21+ GPR34+)", # similar to those identified here https://www.nature.com/articles/s41586-020-2922-4 
               "3" = "Macrophages (HES4+)", # very similar to 5, HES4 is one of its most unique markers
               "4" = "Kupffer cells (SUCNR1+)", # disproves https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1986575/; these KC are inflammatory/antiviral
               "5" = "Macrophages",
               "6" = "Monocytes (TREM2+ CD9+)", # in ramachandran et al, assoc with fibrotic scars
               ## in their projection, also close to KC
               "7" = "cDC2",
               "8" = "cDC1",
               "9" = "Monocytes (secretory)", # related to 2 and 6
               "10" = "activated DCs") # CD80, CD86, CCR7
mk_monocells$annot = new_m_labs[as.character(mk_monocells$cluster)]

imm_df = data.frame(Embeddings(all_imm_cells, reduction = "umap"))
imm_df$all_annot = all_imm_cells$immune_annot
imm_df$simp_annot = imm_df$all_annot
imm_df$simp_annot[imm_df$simp_annot %in% c("ab-T cells (stress)", 
                                           "Hepatocyte-Monocyte interaction")] = "Non annotated cells"
mono = unique(imm_df$all_annot)[c(2,4:9,11:12,31)]
imm_df$simp_annot[imm_df$simp_annot %in% mono] = "Monocytes/Macrophages/cDCs"
tnk = unique(imm_df$all_annot)[c(17:20,22:30,32)]
imm_df$simp_annot[imm_df$simp_annot %in% tnk] = "T cells/ILCs"

cols = MetBrewer::met.brewer("Egypt", length(unique(imm_df$simp_annot)))
names(cols) = unique(imm_df$simp_annot)
cols["Non annotated cells"] = "grey80"

imm_df = imm_df[sample(1:nrow(imm_df), nrow(imm_df), replace = F),]
ggplot()+
  geom_point(aes(x = UMAP_1, y = UMAP_2, colour = simp_annot), 
             imm_df, size = 0.35)+
  guides(colour = guide_legend(override.aes = list(size = 3)))+
  scale_colour_manual(values = cols[order(names(cols))])+
  theme_void()+
  theme(aspect.ratio = 1,
        legend.text = element_text(size = 6),
        legend.title = element_blank())

t_df = data.frame(Embeddings(all_t_cells, reduction = "umap"))
t_df$all_annot = all_t_cells$t_annot

cols = MetBrewer::met.brewer("Juarez", length(unique(t_df$all_annot)))
names(cols) = unique(t_df$all_annot)

t_df = t_df[sample(1:nrow(t_df), nrow(t_df), replace = F),]
ggplot()+
  geom_point(aes(x = UMAP_1, y = UMAP_2, colour = all_annot), 
             t_df, size = 0.35)+
  guides(colour = guide_legend(override.aes = list(size = 3)))+
  scale_colour_manual(values = cols[order(names(cols))])+
  theme_void()+
  theme(aspect.ratio = 1,
        legend.text = element_text(size = 6),
        legend.title = element_blank())

Save data for UMAPs

m_df = data.frame(Embeddings(all_mon_cells, reduction = "umap"))
m_df$all_annot = all_mon_cells$mono_annot

cols = MetBrewer::met.brewer("Klimt", length(unique(m_df$all_annot)))
names(cols) = unique(m_df$all_annot)

m_df = m_df[sample(1:nrow(m_df), nrow(m_df), replace = F),]
ggplot()+
  geom_point(aes(x = UMAP_1, y = UMAP_2, colour = all_annot), 
             m_df, size = 0.35)+
  guides(colour = guide_legend(override.aes = list(size = 3)))+
  scale_colour_manual(values = cols[order(names(cols))])+
  theme_void()+
  theme(aspect.ratio = 1,
        legend.text = element_text(size = 6),
        legend.title = element_blank())

Prepare marker heatmaps

saveRDS(m_df, file = "results/cirrhosis/umap_m_df.RDS")
saveRDS(t_df, file = "results/cirrhosis/umap_t_df.RDS")
saveRDS(imm_df, file = "results/cirrhosis/umap_imm_df.RDS")

Save matrices for myeloid and lymphoid markers

topmk_l = mk_lcells %>%
  group_by(annot) %>%
  top_n(n = 100, wt = avg_log2FC)

Registered S3 method overwritten by 'cli':
  method     from         
  print.boxx spatstat.geom

topmk_t = mk_tcells %>%
  group_by(annot) %>%
  top_n(n = 100, wt = avg_log2FC)
topmk_m = mk_mcells%>%
  group_by(annot) %>%
  top_n(n = 100, wt = avg_log2FC)
topmk_mono = mk_monocells%>%
  group_by(annot) %>%
  top_n(n = 100, wt = avg_log2FC)

mk_use = c("PTGDS","CX3CR1", "KIT", "AHR", "RORC", "CD3E", "CCR7", "MKI67", "CD44", "IL7R",
           "CD80", "CD86", "TREM2", "CD9", "HES4", "IGSF21", "GPR34", "IGHG1", "CD8A", "CD5L",
           "CD4", "FOXP3", "CTLA4", "CD79A", "IGHA1", "IGHM", "CXCL10", "LYVE1", "MARCO",
           "S100A8","LYZ","SELL","IFI27","FCER1A", "ITGA1", "THEMIS","S100B",
           "TRAC", "TRBC1", "TRDC", "TRGC1", "KLRC2", "CCL19", 
           "IDO1", "CLEC9A", "IL2RB", "CXCR6", "CCR6", "CCR5", "RORC", "ZBTB16")

mk_m = mk_use[mk_use %in% topmk_mono$gene | mk_use %in% topmk_m$gene]
mk_l = mk_use[mk_use %in% topmk_l$gene | mk_use %in% topmk_t$gene]

topmk_l = mk_lcells %>%
  group_by(annot) %>%
  top_n(n = 3, wt = avg_log2FC)
topmk_l_l = tapply(topmk_l$gene, topmk_l$annot, function(x) x[1:3])
topmk_t = mk_tcells %>%
  group_by(annot) %>%
  top_n(n = 3, wt = avg_log2FC)
topmk_t_l = tapply(topmk_t$gene, topmk_t$annot, function(x) x[1:3])
topmk_m = mk_mcells%>%
  group_by(annot) %>%
  top_n(n = 3, wt = avg_log2FC)
topmk_m_l = tapply(topmk_m$gene, topmk_m$annot, function(x) x[1:3])
topmk_mono = mk_monocells%>%
  group_by(annot) %>%
  top_n(n = 3, wt = avg_log2FC)
topmk_mono_l = tapply(topmk_mono$gene, topmk_mono$annot, function(x) x[1:3])

rem_ct = c("ab-T cells (stress)", "CD8 ab-T cells (stress)")
id_m = unique(names(c(topmk_m_l, topmk_mono_l)))
id_m = id_m[id_m %in% all_imm_cells$immune_annot & !id_m %in% rem_ct]
feat_m = unique(c(unlist(c(topmk_m_l, topmk_mono_l)[id_m]), mk_m))
feat_m = feat_m[!grepl(".",feat_m,fixed = T)]

id_l = unique(names(c(topmk_l_l, topmk_t_l)))
id_l = id_l[id_l %in% all_imm_cells$immune_annot & !id_l %in% rem_ct]
feat_l = unique(c(unlist(c(topmk_l_l, topmk_t_l)[id_l]), mk_l))
feat_l = feat_l[!grepl(".",feat_l,fixed = T)]

avg_m = AverageExpression(all_imm_cells, features = feat_m, 
                          group.by = "immune_annot")$SCT[,id_m]

Warning in PseudobulkExpression(object = object, pb.method = "average",  :
  Exponentiation yielded infinite values. `data` may not be log-normed.

avg_l = AverageExpression(all_imm_cells, features = feat_l, 
                          group.by = "immune_annot")$SCT[,id_l]

Warning in PseudobulkExpression(object = object, pb.method = "average",  :
  Exponentiation yielded infinite values. `data` may not be log-normed.

mat_m = scale(t(avg_m))
mat_l = scale(t(avg_l))

cols_pal = colorRampPalette(rev(brewer.pal(n = 9, name = "RdBu")))(100)
pheatmap::pheatmap(mat_m, clustering_method = "ward.D", color = cols_pal,
                   treeheight_col = 0, treeheight_row = 20, fontsize = 6.5)

pheatmap::pheatmap(mat_l, clustering_method = "ward.D", color = cols_pal,
                   treeheight_col = 0, treeheight_row = 20, fontsize = 6.5)

Changes in proportions between conditions

saveRDS(mat_m, file = "results/cirrhosis/mat_myeloid_markers.RDS")
saveRDS(mat_l, file = "results/cirrhosis/mat_lymphoid_markers.RDS")

Save proportions data

don_cond_df = unique(all_imm_cells@meta.data[,c("Name", "Donor", "Condition")])
tab_ctCond = table(all_imm_cells$immune_annot, all_imm_cells$Name)
tab_ctCond = tab_ctCond[!grepl("stress", rownames(tab_ctCond)) & 
                          !grepl("intera", rownames(tab_ctCond)),]
imm_props = t(t(tab_ctCond)/colSums(tab_ctCond))*100
imm_props = merge(data.frame(imm_props), don_cond_df, by.x = 2, by.y = 1)

ggplot(imm_props, aes(x = Condition, y = Freq, group = Condition, colour = Condition))+
  facet_wrap(~Var1, scales = "free")+
  geom_jitter(position = position_jitterdodge(jitter.width = 0.3, dodge.width = 1))+
  stat_summary(fun.data = mean_se, position = position_dodge(width = 1), 
               alpha = 0.35, colour = "black")+
  theme_bw()+
  theme(legend.position = "none")

tab_ctCond = table(all_imm_cells$immune_annot, all_imm_cells$Name)
tab_ctCond = tab_ctCond[!grepl("stress", rownames(tab_ctCond)) & 
                          !grepl("intera", rownames(tab_ctCond)),]
imm_cnts = tab_ctCond
imm_cnts = merge(data.frame(imm_cnts), don_cond_df, by.x = 2, by.y = 1)
pvals_l = list()
for(ct in unique(imm_cnts$Var1)){
  dfct = imm_cnts[imm_cnts$Var1==ct,]
  dfct$not = tapply(imm_cnts$Freq, imm_cnts$Var2, sum)-dfct$Freq
  dfct$tot = tapply(imm_cnts$Freq, imm_cnts$Var2, sum)
  dfct$Condition = factor(dfct$Condition, levels = c("healthy", "embolised", "regenerating"))
  mod = glm(cbind(Freq, not) ~ Condition + Donor + tot, data = dfct, family = "binomial")
  pvals_l[[ct]] = summary(mod)$coefficients[c("Conditionembolised", "Conditionregenerating"),
                                            "Pr(>|z|)"]
}

Warning: glm.fit: algorithm did not converge
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

pvals_df = data.frame(t(data.frame(pvals_l)))
rownames(pvals_df) = names(pvals_l)
colnames(pvals_df) = gsub("Condition", "", colnames(pvals_df))
pvals_df$embolised = p.adjust(pvals_df$embolised, method = "fdr")
pvals_df$regenerating = p.adjust(pvals_df$regenerating, method = "fdr")

Load Ramachandran data

saveRDS(imm_props, file = "results/cirrhosis/imm_props_dat.RDS")
saveRDS(pvals_df, file = "results/cirrhosis/imm_props_pval.RDS")

Subset Ramachandran and get HVG

f = "results/cirrhosis/cirr_data_renorm.RDS"
if(!file.exists(f)){
  load("../../data/published/Ramachandran_liver/tissue.rdata")
  
  cirr_data = CreateSeuratObject(counts = tissue@raw.data, meta.data = tissue@meta.data)
  
  # Renormalise, since the original data doesn't look ok
  cirr_data = suppressWarnings(SCTransform(cirr_data, do.correct.umi = T, verbose = F, 
                                           vars.to.regress=c("dataset", "nCount_RNA"),
                                           variable.features.rv.th = 1, seed.use = 1,
                                           return.only.var.genes = F, 
                                           variable.features.n = NULL))
  
  
  cirr_data$anno_cond = paste0(cirr_data$annotation_indepth, "_", cirr_data$condition)
  saveRDS(cirr_data, file = f)
} else{
  cirr_data = readRDS(f)
}

Correlation with Ramachandran data

f = "results/cirrhosis/cirr_data_allimmune_renorm.RDS"
if(!file.exists(f)){
  cirr_imm_data = cirr_data[,cirr_data$annotation_lineage %in% c("MPs", "Tcells", "ILCs",
                                                                 "Bcells", "pDCs", 
                                                                 "Plasma Bcells", "Mast cells")]
  cirr_imm_data = suppressWarnings(SCTransform(cirr_imm_data, do.correct.umi = T, verbose = F, 
                         vars.to.regress=c("dataset", "nCount_RNA"),
                         variable.features.rv.th = 1, seed.use = 1,
                         return.only.var.genes = F, variable.features.n = NULL))
  saveRDS(cirr_imm_data, file = "results/cirrhosis/cirr_data_allimmune_renorm.RDS")
} else{
  cirr_imm_data = readRDS(f)
}

# MPs only
f = "results/cirrhosis/cirr_data_MPs_renorm.RDS"
if(!file.exists(f)){
  cirr_mps_data = cirr_data[,cirr_data$annotation_lineage %in% c("MPs")]
  cirr_mps_data = suppressWarnings(SCTransform(cirr_mps_data, do.correct.umi = T, verbose = F, 
                         vars.to.regress=c("dataset", "nCount_RNA"),
                         variable.features.rv.th = 1, seed.use = 1,
                         return.only.var.genes = F, variable.features.n = NULL))
  saveRDS(cirr_mps_data, file = "results/cirrhosis/cirr_data_MPs_renorm.RDS")
} else{
  cirr_mps_data = readRDS(f)
}

# Tcells and ILCs only
f = "results/cirrhosis/cirr_data_tcells_renorm.RDS"
if(!file.exists(f)){
  cirr_tce_data = cirr_data[,cirr_data$annotation_lineage %in% c("Tcells", "ILCs")]
  cirr_tce_data = suppressWarnings(SCTransform(cirr_tce_data, do.correct.umi = T, verbose = F, 
                         vars.to.regress=c("dataset", "nCount_RNA"),
                         variable.features.rv.th = 1, seed.use = 1,
                         return.only.var.genes = F, variable.features.n = NULL))
  saveRDS(cirr_tce_data, file = f)
} else{
  cirr_tce_data = readRDS(f)
}